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FORMULATION OF THE STATIONARY QUASI-GEOSTROPHIC EQUATIONS 

OF THE OCEAN 



ERICH L FOSTER , TRAIAN ILIESCU , AND ZHU WANG * 

Abstract. This paper presents a conforming finite element discretization of the streamfunction formulation of 
the one-layer stationary quasi-geostrophic equations, which are a commonly used model for the large scale wind- 
driven ocean circulation. Optimal error estimates for this finite element discretization with the Argyris element 
are derived. Numerical tests for the finite element discretization of the quasi-geostrophic equations and two of 
its standard simplifications (the linear Stommel model and the linear Stommel-Munk model) are carried out. By 
benchmarking the numerical results against those in the published literature, we conclude that our finite element 
discretization is accurate. Furthermore, the numerical results have the same convergence rates as those predicted by 
the theoretical error estimates. 
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1. Introduction. With the continuous increase in computational power, complex mathemat- 
ical models are becoming more and more popular in the numerical simulation of oceanic and at- 
mospheric flows. For some geophysical flows in which computational efficiency is of paramount im- 
portance, however, simplified mathematical models are central. For example, the quasi-geostrophic 
equations (QGE), a standard mathematical model for large scale oceanic and atmospheric flows 
[SI HSl [311 [33] , are often used in climate modeling |5] . 

The QGE are usually discretized in space by using the finite difference method (FDM) |32j . 
The finite element method (FEM), however, offers several advantages over the popular FDM, as 
outlined in |30j : (i) an easy treatment of complex boundaries, such as those of continents for the 
ocean, or mountains for the atmosphere; (ii) an easy grid refinement to achieve a high resolution in 
regions of interest [3]; (iii) a natural treatment of boundary conditions; and (iv) a straightforward 
approach for the treatment of multiply connected domains [3D]. Despite these advantages, there 
are relatively few papers that consider the FEM appHed to the QGE [4] [T3 1 [24 l I30l 133] . 

To our knowledge, all the finite element (FE) discretizations of the QGE have been devel- 
oped for the streamfunction-vorticity formulation, none using the streamfunction formulation. The 
reason is simple: The streamfunction-vorticity formulation yields a second order partial differential 
equation (PDF), whereas the streamfunction formulation yields a fourth order PDF. Thus, although 
the streamfunction-vorticity formulation has two variables (g and 4') s-nd the streamfunction for- 
mulation has just one (V'), the former is the preferred formulation used in practical computations, 
since its conforming FE discretization requires low-order (C") elements, whereas the latter requires 
high-order (C^) elements. 

Although the FE discretizations of the QGE are relatively scarce, the corresponding error anal- 
ysis seems to be even more scarce. To our knowledge, all the error analysis for the FE discretization 
of the QGE has been done for the streamfunction-vorticity formulation, and none has been done 
for the streamfunction formulation. Furthermore, to the best of our knowledge, all the available 
error estimates for the FE discretization of the QGE are suhoptimal. The first error analysis for the 
FE discretization of the QGE was carried out by Fix [T3], in which suboptimal error estimates for 
the streamfunction-vorticity formulation were proved. Indeed, relationships (4.7) and (4.8) (and 
the discussion above these) in [T3] show that the FE approximations for both the potential vorticity 
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(denoted by C) and streamfunction (denoted by ip) consist of piecewise polynomials of degree k — 1. 
At the top of page 381, the author concludes that the error analysis yields the following estimates: 



U~i:%^0{h'^-'), (1.1) 

\\C-C% = 0ih'^-'). (1.2) 

Although the streamfunction error estimate ( |1.1[ ) appears to be optimal, the potential vorticity 



error estimate (1.2) is clearly suboptimal. Indeed, using piecewise polynomials of degree k—1 for 
the FE approximation of the vorticity, one would expect an 0{h^) error estimate in the L? norm. 
Medjo [551 [22] used a FE discretization of the streamfunction- vorticity formulation and proved error 
estimates for the time discretization, but no error estimates for the spatial discretization. Finally, 
Gascon et al. [4] proved both a priori and a posteriori error estimates for the FE discretization of 



the linear Stommel-Munk model (see Section 5.2 for more details). This model, while similar to the 



QGE, has one significant difference: the linear Stommel-Munk model is linear, whereas the QGE 
are nonlinear. 

We note that the state-of-the-art in the FE error analysis for the QGE seems to reflect that for 
the two-dimensional Navier-Stokes equations (2D NSE), to which the QGE are similar in form. In- 
deed, as carefully discussed in [1^ (see also [HI [El 1201121] ) , the 2D NSE in streamfunction- vorticity 
formulation are easy to implement (only C° elements are needed for a conforming discretization), 
but the available error estimates are suboptimal (see Section 11.6 in [19]). Next, we summarize the 
discussion in [T^], since we believe it sheds light on the QGE setting. For C° piecewise polynomial 
of degree k FE approximation for both the vorticity (denoted by w) and streamfunction (denoted 
by ip), the error estimates given in fTT] are (see (11.26) in [19]): 

IV- - V'^li + l|w - ^''llo < Ch''-^/^ I In/ir, (1.3) 

where a ~ 1 for k — 1 and cr for fc > 1. It is noted in [TS] that the error estimate in ( |1.3| ) is 
not optimal: one may loose a half power in h for the derivatives of the streamfunction (i.e., for the 
velocity), and three-halves power for the vorticity. It is also noted that there is computational and 



theoretical evidence that (1.3) is not sharp with respect to the streamfunction error. Furthermore, 
in |14j it was shown that, for the linear Stokes equations, the derivatives of the streamfunction are 
essentially optimally approximated (see (11.27) in [W): 

l^.^/^l^ <C/^^-^ (1.4) 



where e — for k > 1 and e > is arbitrary for fc — 1. It is, however, noted in [19] that ( [1.3[ ) seems 
to be sharp for the vorticity error and thus vorticity approximations are generally poor. 

The streamfunction formulation is, from both mathematical and computational points of view, 
completely different from the streamfunction-vorticity formulation. Indeed, the FE discretization 
of the streamfunction formulation generally requires the use of elements (for a conforming 
discretization), which makes their implementation challenging. From a mathematical point of 
view, however, the streamfunction formulation has the following significant advantage over the 
streamfunction-vorticity formulation: there are optimal error estimates for the FE discretization of 
the streamfunction formulation (see the error estimate (13.5) and Table 13.1 in [IS]), whereas the 
available error estimates for the streamfunction-vorticity formulation are suboptimal. 

The main goal of this paper is twofold. First, we use a finite element (the Argyris element) to 
discretize the streamfunction formulation of the QGE. To the best of our knowledge, this is the first 
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time that a finite element has been used in the numerical discretization of the QGE. Second, 
we derive optimal error estimates for the FE discretization of the QGE and present supporting 
numerical experiments. To the best of our knowledge, this is the first time that optimal error 
estimates for the QGE have been derived. 

The rest of the paper is organized as follows: Section [2] presents the QGE, their weak formula- 
tion, and mathematical support for the weak formulation. Section[3]outlines the FEM discretization 
of the QGE, posing a special emphasis of the Argyris element. Rigorous error estimates for the 
FE discretization of the stationary QGE are derived in Section |4j Several numerical experiments 
supporting the theoretical results are presented in Section [5j Finally, conclusions and our future 
research directions are included in Section [6l 

2. The Quasi- Geostrophic Equations. The large scale ocean flows, which play a significant 
role in climate dynamics [21 [H], are driven by two major sources: the wind and the buoyancy (see, 
e.g.. Chapters 14-16 in Winds drive the subtropical and subpolar gyres, which correspond to 

the strong, persistent, subtropical and subpolar western boundary currents in the North Atlantic 
Ocean (the Gulf Stream and the Labrador Current) and North Pacific Ocean (the Kuroshio and 
the Oyashio Currents), as well as their subtropical counterparts in the southern hemisphere [^135). 
One of the common features of these gyres is that they display strong western boundary currents, 
weak interior flows, and weak eastern boundary currents. 

One of the most popular mathematical models used in the study of large scale wind-driven 
ocean circulation is the QGE [H [35] . The QGE represent a simplified model of the full-fledged 
equations (e.g., the Boussinesq equations), which allows efficient numerical simulations while pre- 
serving many of the essential features of the underlying large scale ocean flows. The assumptions 
used in the derivation of the QGE include the hydrostatic balance, the /3-plane approximation, the 
geostrophic balance, and the eddy viscosity parametrization. Details of the derivation of the QGE 
and the approximations used along the way can be found in standard textbooks on geophysical 
fluid dynamics, such as [H [25l [Ml [23 [SH [35] . 

In the one-layer QGE, sometimes called the barotropic vorticity equation, the flow is assumed to 
be homogenous in the vertical direction. Thus, stratification effects are ignored in this model. The 
practical advantages of such a choice are obvious: the computations are two-dimensional, and, thus, 
the corresponding numerical simulation have a low computational cost. To include stratification 
effects, QGE models of increasing complexity have been devised by increasing the number of layers 
in the model (e.g., the two-layer QGE and the A''-layer QGE [3S]). As a first step, in this report 
we use the one-layer QGE (referred to as "the QGE" in what follows) to study the wind-driven 
circulation in an enclosed, midlatitude rectangular basin, which is a standard problem, studied 
extensively by ocean modelers [H [23 IMl (UliSIl (33 • 

The nondimensional streamfunction-vorticity formulation of the stationary one-layer quasi- 
geostrophic equations is (see, e.g., equation (14.57) in [35] , equation (1.1) in |26| . equation (1.1) in 
[5S] . and equation (1) in [T5]): 

J{^,q) = -Re-^Aq + F (2.1) 
q = -RoAtJj + y, (2.2) 

where ^ is the velocity streamfunction, q is the potential vorticity, F is the forcing, J(-,-) is the 
Jacobian operator given by 

d^p dq dip dq 



Re is the Reynolds number, and Ro is the Rossby number. The Rossby number, i?o, is defined as 



Ro 



U 



(2.4) 



where /3 is the coefhcient multiplying the y coordinate in the /3-plane approximation [S] [3S] , L is 
the width of the computational domain, and U is the Sverdrup velocity obtained from the balance 
between the /3-effect and the curl of the divergence of the wind stress [3S1. The Reynolds number, 
Re, is defined as 



Re := 



UL 



(2.5) 



where A is the eddy viscosity parametrization. The horizontal velocity u can be recovered from tp 

( dip dtp 
dy ' dx 



by using the following formula: u = ( ^ 



Substituting (2.2) in (2.1) and dividing by Ro, we get the streamfunction formulation of the 



stationary one-layer quasi-geostrophic equations 

+ J(V',AV') 



Re- 



Ro-^^^Ro-'F. 
ox 



(2.6) 



We note that the streamfunction- vorticity formulation has two unknowns {q and whereas 
the streamfunction formulation has only one unknown (ip). Because the streamfunction- vorticity 
formulation is a second-order PDE, whereas the streamfunction formulation is a fourth-order PDE, 
the former is more popular in practical computations. 

We also note that (2.1)-(2.2| and (2.6) are similar in form to the 2D NSE written in the 
streamfunction-vorticity and streamfunction formulations, respectively. There are, however, sev- 



eral significant differences between the QGE and the 2D NSE. First, the term y in (2.2) and the 
corresponding term ^ in (2.6), which model the rotation effects in the QGE, do not have coun- 
terparts in the 2D NSE. Second, the Rossby number, Ro, in the QGE, which is a measure of the 
rotation effects, does not appear in the 2D NSE. 

Next, we comment on the significance of the two parameters in (2.6), the Reynolds number. 
Re, and the Rossby number, Ro. As in the 2D NSE case. Re is the coefficient of the diffusion term 
—Aq = A'^ip. The higher the Reynolds number Re, the smaller the magnitude of the diffusion term 
as compared with the nonlinear convective term Ji^ip, Aip). For small Ro, which corresponds to large 
rotation effects, the forcing term, Ro^^ F, becomes large compared with the other terms. The term 



Ro- 



-1 d± 

dx 



could be interpreted as a convection type term with respect to ip, not to g = —Ai/j. When 



Ro is small, Ro~^ ^ becomes large. Thus, the physically relevant cases for large scale oceanic flows, 
in which Re is large and Ro is small (i.e., small diffusion and high rotation, respectively) translate 
mathematically into a convection-dominated PDE with large forcing. Thus, from a mathematical 
point of view, we expect the restrictive conditions used to prove the well-posedness of the 2D NSE 
[Tni[T71[TH] to be even more restrictive in the QGE setting, due to the rotation effects. We will later 
see that this is indeed the case. 

To completely specify the equations in (2.6), we need to impose boundary conditions. The 
question of appropriate boundary conditions for the QGE is a thorny one, especially for the 
streamfunction-vorticity formulation (see, e.g., OEl]). In this report, we consider ip — ^ = on 
dft, which are also used in 021 for the streamfunction formulation of the 2D NSE. 



To derive the weak formulation of the QGE (2.6), we first introduce the appropriate functional 



setting. Let X := H^{^) = jV' G H'^i^) : -0 = |^ = on (9f}|. Multiplying by a test function 
X X and using the divergence theorem, we get the weak formulation of the QGE in strcamfunction 
formulation 1191: 



i^xXdx 



Therefore, letting 



ao(V',x) 


= Re-' 




-h 




Jn 


a2{ip,x) 


= -Ro 


eix) 


= Ro-' 



(2.7) 

(2.8) 
(2.9) 
(2.10) 
(2.11) 

gives the weak formulation of the QGE in strcamfunction formulation: Find ijj (z X such that 

ao(V',x)+ai(^,^,x)+a2(V',x) = ^(x), Vx G X (2.12) 

The linear form £, the bilinear forms ao and 02, and the trilinear form ai are continuous: There 
exist Fi > and F2 > such that 



Alp Ax dx, 
AC {ipyXx -ipxXy) dx, 
X dx, 
Fxdx, 



|ao(^,x)l <i?e-i|^l2|xl2 V^,xG^, 

|ai(C,V',x)l <ri|Cl2|V'|2|xl2 vc,0,x€X, 

|a2(^,x)l <i?o-iF2|V'|2|x|2 V0, xe^, 

K(x)| <i?o"' 11^11-2 |x|2 Vxex 



(2.13) 
(2.14) 
(2.15) 
(2.16) 



Inequalities (2.13), (2.14), and (2.16) are stated in |5] (see inequalities (2.2) and (2.3) in [5") 
Inequality (2.15) can be proved as follows. Proposition 2.1(iii) in [28J implies that 



|a2(V',x)l <i?o"'C 11^112 ||xl|2, 
where C is a generic constant. Theorem 1.1 in [TT implies that 



(2.17) 



2, the iJ^-seminorm, and 



the iJ^-norm are equivalent on X = Hq. Thus, (2.17) yields inequality ( 2.15| ). 

For small enough data, one can use the same type of arguments as in [161 [T7] to prove that the 
QGE in streamfunction formulation (2.12) are well-posed [J 137]. In what follows, we will always 
assume that the small data condition involving Re, Ro and F, is satisfied and, thus, that there 



exists a unique solution -0 to (2.12) 



Using a standard argument [5 , one can also prove the following stability estimate: 



Theorem 2.1. The solution ip of (2.12) satisfies the following stability estimate: 

I0I2 < ReRo-' ||F||_2. 



(2.18) 



Proof. Setting x = -0 in (2.12 ), wc get: 



ao(V',V')+ai(^,^,V') + a2(^,^) ^^W- (2.19) 
Since the trilinear form ai is skew-symmetric in the last two arguments |16L I17L 119) , we have 

ai(V^»)=0. (2.20) 
We also note that, applying Green's theorem, we have 



Rq-^ 



dx 



ip dx dy 



Rq- 



d_ 
dx 



(■0^) dx dy 



d_ 
dx 



dy 



(0) dx dy 



Ro~ 



Odx + Tp'^dy = 0, (2.21) 



on 



where in the last equality in (2.21) we used that -0 = on dQ (since -0 G Hq{Q)). Substituting 
(2.21) and (2.20) in (2.19) and using the Cauchy-Schwarz inequality, we get: 

Ai: Alp dyi ^ Re Rq-^ / F ip d:ic < Re Ro^^ \\F\\^2\'4'\2, 
Jn 



ml 



(2.22) 



which proves (2.18). 



□ 



3. Finite Element Formulation. In this section, we present the functional setting and some 
auxiliary results for the FE discretization of the streamfunction formulation of the QGE (2.12). Let 
T'* denote a finite element triangulation of ft with meshsize (maximum triangle diameter) h. We 
consider a conforming FE discretization of (2.12), i.e., X''' C X = Hq{Q). 



The FE discretization of the streamfunction formulation of the QGE ( 2.12 ) reads: Find 0'* € X'^ 
such that 

ao(0^x'^)+al(V'^V'^x'^)+a2(V'^x") = ^(x"), vx^ex". (s.i) 

Using standard arguments |161 117j. one can prove that, if the small data condition used in proving 



the well-posedness result for the continuous case holds, then (3.1) has a unique solution ■0'* (see 
Theorem 2.1 and subsequent discussion in [5J). One can also prove t he fo llowing stability result for 
'0'' using the same arguments as those used in the proof of Theorem 2.1 for the continuous setting. 



Theorem 3.1. The solution of (3.1) satisfies the following stability estimate: 



\iP%<ReRo-^\\F\\ 



(3.2) 



In order to develop a conforming FEM for the QGE (2.12), we need to construct subspaces of 
Hq(0), i.e., to find EEs, such as the Argyris triangular element, the Bell triangular element, the 
Hsieh-Clough- Tocher triangular element (a macroelement), or the Bogner-Fox-Schmit rectangular 
element [1311213]. In what follows, we will use the Argyris FE. The Argyris FE employs piecewise 
polynomials of degree five and has twenty-one degrees of freedom (DOFs): the value at each vertex, 
the value of the first derivatives at each vertex, the value of the second derivatives at each vertex. 



the value of the mixed derivative at each vertex, and the value of the normal derivatives at each 
of the edge midpoints. To maintain the direction of the normal derivatives in the transformation 
from the reference element to the physical element, we use the approach developed in |10j . 

By using Theorem 6.1.1 and inequality (6.1.5) in we obtain the following three approxima- 
tion properties for the Argyris FE space X'^: 



yx(^ H\n)nH'^{n), ax'^ex'' 



such that 
such that 
such that 



\\x~-X%<Ch^\x\e. 

\\x-x%<Ch^\xU, 
Ilx-x"ll2 <c;i|x|3, 



(3.3) 
(3.4) 
(3.5) 



where C is a generic constant that can depend on the data, but not on the meshsize h. Property 
(3.3) follows from (6.1.5) in |S] with q = 2, p — 2, m 
(6.1.5) in [Sj with g = 2, p = 2, to = 2 and fc + 1 = 
in [6 with q = 2, p — 2, 



6. Property (3.4) follows from 



= 2 and fc + 1 
4. Finally, property (3.5) follows from (6.1.5) 



2 and fc + 1 



4. Error Analysis. The main goal of this section is to develop a rigorous numerical analysis 



In 



for the FE discretization of the QGE (3.1 ) by using the conforming Argyris element. In Theorem |4.1 
we prove error estimates in the norm by using an approach similar to that used in [S]. 
we prove error estimates in the and 
Theorem 4.1 



Theorem 



norms by using a duality argument. 
Let ■0 be the solution of (2.12) and -0'' 



be the solution of (3.1). Furthermore, 



assume that the following small data condition is satisfied: 



Re^^ Ro > Fi 



\F\\. 



(4.1) 



where Re is the Reynolds number defined in (2.5), Ro is the Rossby number defined in (2.4), F 



is the continuity constant of the trilinear form oi in (2.14), and F is the forcing term. Then the 
following error estimate holds: 



IV'-V'^b < C(ii'e,i?o,Fi,F2,F) inf \i, - 



h\ 

X b, 



where T2 is the continuity constant of the bilinear form 02 in (2.15) an 

Ro-^r 



C(i?e,i?o, Fi,F2,F) := 



2 ^2Re-' + ri Re Rq-' ||F||_2 
i?e-i - TiReRo-^ IIFII 2 



(4.2) 



(4.3) 



is a generic constant that can depend on Re, Ro, Fi, F2, F, but not on the meshsize h. 

Remark 4.2. Note that the small data condition in Theorem \4.1\ involves both the Reynolds 
number and the Rossby number, the latter quantifying the rotation effects in the QGE. 

Furthermore, note that the standard small data condition Re^^ > Fi ||i^||_2 used to prove the 
uniqueness for the steady-state 2D NSE il6i. \23l \34^ is significantly more restrictive for the 
QGE, since (4.1) has the Rossby number (which is small when rotation effects are significant) on 
the left-hand side. This is somewhat counterintuitive, since in general rotation effects are expected 
to help in proving the well-posedness of the system. We think that the explanation is the following: 
Rotation effects do make the mathematical analysis of 3D flows more amenable by giving them a 
2D character. We, however, are concerned with 2D flows (the QGE). In this case, the small data 



condition (4.1) (needed in proving the uniqueness of the solution) indicates that rotation effects 



make the mathematical analysis of the (2D) QGE more complicated than that of the 2D NSE. 
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Proof. Since X'' C X, ( |2.12[ ) holds for all x = x'' G X''. Subtracting ([sl]) from ( |2.12[ ) with 
X = x'' G -'^'^ gives 

ao(V - x') + ai(VA V-, X'^) - x") + a^O^ - x') =0 Vx'' e X''. (4.4) 



Next, adding and subtracting ai(-(/''', V'l x'') to (4.4), we get 

ao(^ - ^\x'') + ai(V',V',x") - al(V'^ V,X'') +«l(^^V',x'') " ai(7^\ V^x") 

+a2(^-V'',x') = Vx''eX^ (4.5) 

The error e can be decomposed as e := ip — 4''^ = {ip — A'*) + (A'^ — V''') ?? + ^fi^, where A'' G X'' 
is arbitrary. Thus, equation (4.51 can be rewritten as 

aoiv + V^ x'') + ai{v + ^^ V, x") + + ^^ x') + «2(r7 + ^^ X'^) = O Vx" e X\ (4.6) 

Letting x'' '■— in ( |4.6[ ), we obtain 

ao(v5^ ¥>") + a2(^'\ ^") - -00(^7, ¥>") - ai(??, ^, ¥>") - ai(^">, 



(4.7) 



Note that, since a-2.y^,Lp^) ^ -a2{(p^' , ifi'') V(^'' e X^' C X = H^, it follows that a2((^'', (ys'*) = 0. 
We also have that ai{ip'^, (p'^, (^'') = 0. Using these equalities in (4.7), we get 

ao(^^ ^") = -aoiv, v"") - «i(r?, ^, v"") - a,{^\ij, ^'') - al(V'^ 77, ^") - 02(77, p''). (4.8) 
Using aQ{ip'^,ip''^) — Re^^ W'^W ^^'^ ( |2.8[ ) - ( |2.10[ ) in (4.8), simplifying, and rearranging terms, gives 
\ip% < {Re'^ - Ti I7AI2)"' (i?e-i + Ti |7/;|2 + Ti + Ro'^ T2) hb- (4.9) 
Using (4.91 and the triangle inequality along with the stability estimates ( 2.18[ ) and (3.2 1, gives: 



|e|2 < |?7|2 + I'^'^b < 



i?e-i + Ti IVI2 + Ti + i?o-i r2 



iio-iPa + 2 i?e-i +riEe i?o-i ||F||_2 
Re-^ - TiReRo-^ \\F\\^2 



(4.10) 



where A'' S X'' is arbitrary. Taking the infimum over A'' S X'' in (4.10) proves estimate (4.2). □ 
Next, we prove error estimates in the norm and seminorm by using a duality argument. 

(4.11) 



To this end, we first notice that the QGE (2.6) can be written as 

Mtp^Ro^^ F, 



where the nonlinear operator M is defined as 



Nip:^ Re-' A> + Jii^, AiP) - Ro 



dx ' 

The linearization of M around ijj, a solution of ( |2.6[ ), yields the following linear operator: 

-1 dx 



£ X := Re-' A^x + J(x, AtA) + J{^, Ax) - Ro' 
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dx ' 



(4.12) 



(4.13) 



To find the dual operator C* of £, we use (4.131 and apply Green's theorem: 



{Cx,r)^[ Re-' A'x + Jix, AV) + J(^, Ax) - Ro 



-1 dx 
dx 



X , Re-' tP* - J{^, A^P*) + Ro 



-I or 

dx 



r 



+ lx,J{Ai,,r)] =(x,/:*v^*).(4.i4) 



Thus, the dual operator C* is given by 



C* i)* = Re-^ A^ V* - J{i^, AV/) + JiAiP, iP*) + Ro 
For any given g e L^(rt), the weak formulation of the dual problem is: 

{£*r,x) = {9,x) ^xeX = Hl{^). 



dx 



(4.15) 



(4.16) 



We assume that ip* , the solution of (4.16), satisfies the following elliptic regularity estimates: 



\rh<C\\g\\o, 
\rh<C\\g\\- 



(4.17) 
(4.18) 
(4.19) 



where C is a generic constant that can depend on the data, but not on the meshsize h. 

Remark 4.3. We note that this type of elliptic regularity was also assumed in for the 
streamfunction formulation of the 2D NSE. In that report, it was also noted that, for a polygonal 
domain with maximum interior vertex angle < 126°, the assumed elliptic regularity was actually 
proved in f^. We note that the theory developed in 0/ carries over to our case. In Section 5 in 
|2/ it is proved that, for weakly nonlinear problems that involve the biharmonic operator as linear 
main part and that satisfy certain growth restrictions, each weak solution satisfies elliptic regularity 



results of the form (4.17)-(4.19|. Assuming that ft is a bounded polygonal domain with inner angle 

Theorem 7 in JB/ with k — and k — 1 
j/ to prove that 



uj at each boundary corner satisfying oj < 126.283696 . . 

Using an argument similar to that used in Section 6(b) in 



implies (4.17)-(4.19) 



the streamfunction formulation of the 2D NSE satisfies the restrictions in Theorem 7, we can prove 
that ip* , the solution of our dual problem (4.16), satisfies the elliptic regularity results in (4.17)- 



(4.19). Indeed, the main point in Section 6(b) in ^ is that the corner singularities arising in flows 
around sharp corners are essentially determined by the linear main part in the streamfunction 



formulation of the 2D NSE, which is the linear main part of our dual problem (4.16) as well. 

Theorem 4.4. Let ip be the solution of (2.12) and %p^ be the solution of (3.1). Assume that 



the same small data condition as in Theorem \4.I\ is satisfied. 



Re-^Ro > Ti IIFII 2 



(4.20) 



Furthermore, assume that ip € H^{fl) n i?g(r2). Then there exist positive constants Cq, Ci and C2 
that can depend on Re, Ro, Ti, T2, F, but not on the meshsize h, such that 

\-p-i^%<C2h\ (4.21) 
\-p-^\<Cih\ (4.22) 
U-^"h<C^h\ (4.23) 
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Remark 4.5. The Argyris FE error estimates in Theorem 4-4 can be extended to other con- 
forming FE spaces. 

Proof. Estimate (4.211 follows immediately from (3.3) and Theorem 4.1 Estimates (4.23) and 



(4.22) follow from a duality argument. 

The error in the primal problem (2.12) and the interpolation error in the dual problem (4.16) 
(with the function g to be specified later) are denoted as e := ij^—ip'^ and e* := 'ip*—ip*^, respectively. 

To prove the norm estimate (4.23), we consider g — e in the dual problem (4.16): 

|e|2 = (e,e) = (£e,V*) = ie,C*i^*) - (e, £* e*) + (e, £* V*'') = (£ e, e*) + (£ e>*''). (4.24) 



The last term on the right-hand side of (4.24 1 is given by 



(£ e, V ) = Re~' A^e + J(e, A V) + J(?A, A e) - Ro' 



.1 de 
dx 



(4.25) 



To estimate this term, we consider the error equation obtained by subtracting (3.1 ) (with ?/;'' — ^*'^) 
from ([2I2I) (with X = -0*''): 



Re-'A^e-Ro~'^ 
ox 



v^*"] + (j(V',Av)- J(V'^Av''), v^*") 



Using (4.26), equation (4.25) can be written as follows: 

(£ e, r'') = ( J(e, A 7/;) + J(V', A e) - J(V', A V) + J(^'\ A ^'^) , V*'*) • 

Thus, by using ( |4.27 1 equation (4.24) becomes 

\e\' = iCe,e*) + iCe,r'') 

= ao(e, e*) + 02(6, e*) + ai(e,V', e*) + ai(?/>, e,e*) + ai(e, V', V'*'') 

+ai(V, e, V'*'') - ai(^, ^, ^*'') + al(V'^ V*') 
= ao(e, e*) + 02(6, e*) + ai(e,V', e*) + ai{tp, e,e*) 

-ai(e,V',e*) + ai(e, V''', e*) + ai{e,e,ip*) 



(4.26) 



(4.27) 



(4.28) 



Using the bounds in ( |2.13[ )-( |2.15[ ), ( |4.28| ) yields 

|e|2 < i?e-i |e|2 |e*|2 + i?o-i r2 |e|2 |e*|2 + Ti |e|2 IVI2 |e*|2 + Ti IVI2 \e*\2 
+ri |e|2 |^|2 |e*|2 + Ti |e|2 IV'^b |e*|2 + Ti |e|2 |e|2 
= |e|2 |e*|2 (-Re-i + Ro-^ r2 + Ti |i/^|2 + Ti IV'b + Ti |^|2 + Fi [V^''^) + \e\l {T, \r I2) . (4.29) 



Using the stability estimates (2.181 and (3.2|, (4.29) becomes 

\e\'<C\e\2\e*\2 + \e\l {T.irU), (4.30) 

where C is a generic constant that can depend on Re, Ro, Fi, F2, F, but not on the meshsize h. 
Using the approximation results (3.4), we get 



\e*\2<Ch'\r\4. 
10 



(4.31) 



Using (4.17)-(4.18), the elliptic regularity results of the dual problem (4.16) with g :— e, we get 



1^-14 <C|e|, 



(4.32) 



which obviously implies 



Inequalities (4.31 )-(4.32 1 imply 



\r\2<C\e\ 



< Ch'\e\ 



Inserting ( |4.33[ ) and ( |4.34| in ( |4.30[ ), we get 

|ep<C/i2|e|2|e|+C|e|2|e|. 



Using the obvious simplifications and the H error estimate (4.211 in (4.35) yields 



< Ch^ \e\2 + C\e\l < C + C = Coh'^ 



(4.33) 



(4.34) 



(4.35) 



(4.36) 



which proves the L error estimate (4.231 



Estimate (4.22) can be proven using the same duality argument as that used to prove estimate 
(4.23). The major differences are that we use g — — Ae in the dual problem (4.16) and we use the 
approximation result (3.5). □ 



5. Numerical Results. The main goal of this section is twofold. First, we show that the 



FE discretization of the streamfunction formulation of the QGE (3.1) with the Argyris element 



produces accurate numerical approximations, which are close to those in the published literature 
[H [301 ES]. Second, we show that the numerical results follow the theoretical error estimates in 
[Theorem 4.11 and [Theorem 4.41 

5.1. Mathematical Models. Although the pure streamfunction formulation of the steady 



QGE ( 2.6 ) is our main concern, we also test our Argyris FE discretization on two simplified settings: 
(i) the Linear Stommel model; and (ii) the Linear Stommel-Munk model. The reason for using 
these two additional numerical tests is that they are standard test problems in the geophysical fluid 
dynamics literature (see, e.g., Chapter 14 in Vallis [35] as well as the reports of Myers and Weaver 
[50] and Gascon et al. [3]). This allows us to benchmark our numerical results against those in the 
published literature. Since both the Linear Stommel and the Linear Stommel-Munk models lack 



the nonlinearity present in the QGE (2.6), they represent good stepping stones for testing our FE 
discretization. 

The Linear Stommel-Munk model (see equation (14.42) in [35 and Problem 2 in 4 ) is 



(5.1) 



The parameters eg and tM in ( |5.1[ ) are the Stommel number and Munk scale, respectively, which 
are given by (see, e.g., equation (10) in [30] and equations (14.22) and (14.44) in [S^) eM — 
and es = where A is the eddy viscosity parameterization, /? is the coefficient multiplying the y 
coordinate in the /3-plane approximation, L is the width of the computational domain, and 7 is the 
coefficient of the linear drag (Rayleigh friction) as might be generated by a bottom Ekman layer 
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(see equation (14.5) in [35 ). The model is supplemented with appropriate boundary conditions, 
which will be described for each of the subsequent numerical tests. 

We note that the Linear Stommel-Munk model ( 5.1 ) is similar in form to the QGE ( 2.6 ). Indeed, 



both models contain the biharmonic operator A^-)/), the rotation term and the forcing term /. 
The two main differences between the two models are the following: First, the QGE are nonlinear, 
since they contain the Jacobian term J{tp,q), whereas the Stommel-Munk model is linear. The 
second difference is that the Linear Stommel-Munk model contains a Laplacian term A^, whereas 
the QGE do not. 

We also note that the two models use different parameters: the Reynolds number. Re, and the 
Rossby number, Ro, in the QGE and the Stommel number, eg, and the Munk scale, e^, in the 
Linear Stommel-Munk model. The parameters cm, Ro, and Re are related through em = RoRe~^. 
There is, however, no explicit relationship among es, Ro, and Re. The reason is that the QGE 
(2.6) do not contain the Laplacian term that is present in the Stommel-Munk model (5.1), which 
models the bottom Rayleigh friction. Thus, the coefficient 7 does not have a counterpart in the 
QGE. This explains why £5, which depends on 7, cannot be directly expressed as a function of Ro 
and Re. 

The second simplified model used in our numerical investigation is the Linear Stommel model 
(see, e.g., equation (14.22) in [3S] and equation (11) in [3U|): 



dx 



(5.2) 



We note that the Linear Stommel model (5.2) is just the Linear Stommel-Munk model (5.11 in 
which the biharmonic term is dropped (i.e., cm = 0). 



5.2. Numerical Tests. In this section, we present results for the Linear Stommel model (5.2 1, 



the Linear Stommel-Munk model (5.1), and the (nonlinear) QGE (2.6) 



5.2.1. Linear Stommel Model. This section presents the results for the EE discretization 
of the Linear Stommel model (5.2) by using the Argyris element. The computational domain is 
Q. — [0,1] X [0,1]. For completeness, we present results for two numerical tests. The first test, 
denoted by Test 1, corresponds to the exact solution used by Vallis (equation (14.38) in [3S]), while 
the second test, denoted by Test 2, corresponds to the exact solution used by Myers and Weaver 
(equations (15) and (16) in [30]). 

Test la: In this test, we choose the same setting as that used in equation (14.38) in [3S]. In 
particular, the forcing term and the non-homogeneous Dirichlet boundary conditions are chosen to 
match those given by the exact solution il){x,y) = {1 — x — e~^/'s) gin (ttj/). We choose the same 
Stommel number as that used in 1351 . i.e., eg = 0.04. 



Figure 5.1(a) presents the streamlines of the approximate solution obtained by using the Argyris 



element on a mesh with h ^ ^ and 9670 DoFs. We note that Figure 5.1(a) resembles Figure 14.5 



in |35j . Since the exact solution is available, we can compute the errors in various norms. Table 5.1 
presents the errors eo, ei, and 62 (i.e., the L^, , and H 

of the meshsize, h (the DoFs are also included). We note that the errors in Table 5.1 



theoretical rates of convergence predicted by the estimates (4.21 )-(4.23 ) in Theorem 4.4 



errors, respectively) for various values 

follow the 
The orders 



of convergence in Table 5.1 are close to the theoretical ones for the fine meshes, but not as close for 
the coarse meshes. We think that the inaccuracies on the coarse meshes are due to their inability 
to capture the thin boundary layer at x = 0. The finer the mesh gets, the better this boundary 
layer is captured and the better the numerical accuracy becomes. 
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h 


DoFs 


eo 


L2 order 


ei 


order 


£2 


H'' order 


1/2 


70 


0.1148 




1.81 




83.67 




1/4 


206 


0.01018 


3.495 


0.312 


2.537 


25.48 


1.716 


1/8 


694 


0.0004461 


4.512 


0.02585 


3.593 


3.902 


2.707 


l/l6 


2534 


1.09 X 19~5 


5.355 


0.001215 


4.412 


0.3494 


3.481 


1/32 


9670 


1.972 X 19"^ 


5.788 


4.349 X 19"'^ 


4.804 


0.02335 


3.903 



Table 5.1 




Test lb: To verify whether the degrading accuracy of the approximation is indeed due to 
the thin (western) boundary layer, we use eg = 1 in Test la, which will result in a much thicker 
western boundary layer. We then run Test la, but with the new es- As can be seen in Tabic 5.2 
the rates of convergence are the expected theoretical orders of convergence. This shows that the 
reason for the inaccuracies in [Table 5.1| were indeed due to the thin western boundary layer. 

Test 2: For this test, we use the exact solution given by equations (15) and (16) in [50], i.e., 



:i(7ry) 



- x(l+47r2e|) 



|27res sin(7ra;) + cos(7rx) 



[(1 + e^2)e'"i^ - (1 + e'"i)e^2^] | , where 



Ri^2 = ^^^2e^'*^ "'^ ■ The forcing and the homogeneous Dirichlet boundary conditions are chosen 
to match those given by the exact solution. We choose the same Stommel number as that used in 
[30], i.e., es = 0.05. 

[Figure 5.2] presents the streamlines of the approximate solution obt ained by using the Argyris 
element on a mesh with h = 



32 



and 9670 DoFs. We note that Figure 5.2 resembles Figure 2 in 
[30| . Table 5.3 presents the errors eo, ei, and 62 for various meshsizes h. The errors in Table 5.3 



follow the theoretical rates of convergence predicted by the estimates (4.21 1 - (4.23) in Theorem 4.4 



Again, we see that the orders of convergence in [Table 5.3| are close to the theoretical ones for the 
fine meshes, but not as close for the coarse meshes. We again attribute this to the inaccuracies at 
the thin (western) boundary layer at x = 0. 

5.2.2. Linear Stommel-Munk Model. This section presents results for the FE discretiza- 
tion of the Linear Stommel-Munk model ( [5.1[ ) by using the Argyris element. Our computational 
setting is the same as that used by Gascon et al. [4]: The computational domain is = [0, 3] x [0, 1], 
the Munk scale is ejv/ = 6 x 10~^, the Stommel number is eg = 0.05, and the boundary conditions 
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h 


DoFs 


eo 


L2 order 


ei 


order 


£2 


H'^ order 


1/2 


70 


1.689 X 10"^ 




0.0003434 




0.008721 




1/4 


206 


3.722 X 10"^ 


5.504 


1.341 X 10"^ 


4.678 


0.0005616 


3.957 


1/8 


694 


4.891 X 10"^ 


6.25 


3.757 X 10"^ 


5.158 


3.25 X 10"^ 


4.111 


l/l6 


2534 


7.079 X 10^" 


6.111 


1.117 X 10"** 


5.071 


1.964 X 10"" 


4.049 


1/32 


9670 


1.08 X 10~^^ 


6.035 


3.437 X 10"^° 


5.023 


1.213 X 10"^ 


4.018 



Table 5.2 

Linear Stommel Model \5.2\, Test lb \35f : The errors eo, ei, 62 for various meshsizes h. 



h 


DoFs 


eo 


L2 order 


ei 


order 


62 


H'' order 


1/2 


70 


0.005645 




0.1451 




6.602 




1/4 


206 


0.0004276 


3.723 


0.02081 


2.801 


1.632 


2.016 


1/8 


694 


1.46 X lO"'^ 


4.872 


0.001408 


3.886 


0.2066 


2.982 


I/16 


2534 


2.954 X 10"'^ 


5.627 


5.829 X 10-5 


4.594 


0.0165 


3.646 


1/32 


9670 


4.968 X 10-9 


5.894 


1.998 X 10-*^ 


4.867 


0.001069 


3.948 



Table 5.3 

Linear Stommel Model \5.2\, Test 2 \30^: The errors eo, e\, e2 for various meshsizes h. 



are = ^ = on dVt. For completeness, we present results for two numerical tests, denoted by 
Test 3 and Test 4, corresponding to Test 1 and Test 2 in [5, respectively. 

Test 3: For this test, we use the exact solution given by Test 1 in [3], i.e., ip{x,y) = 
sin^ (^) sin^ (ttj/). The forcing term is chosen to match that given by the exact solution. 



Figure 5.3(a)| presents the streamlines of the approximate solution obtained by using the Argyris 



element on a mesh with h = ^ and 28550 DoFs. We note that Figure 5.3(a) resembles Figure 7 in 



30] . Table 5.4] presents the errors eo, ei, and 62 for various meshsizes h. The errors in Table 5.4 



follow the theoretical rates of convergence predicted by the estimates ( 4.21[ )-(4.23) in Theorem 4.4 



This time, we see that the orders of convergence in [Table 5.4| are close to the theoretical ones for 
the fine meshes, but are higher than expected for the coarse meshes. We attribute this to the fact 
that the exact solution does not display any boundary layers that could be challenging to capture 
by the Argyris element on a coarse mesh. 

Test 4: For this test, we use the exact solution given by Test 2 in [3], i.e., ip{x,y) = 
[(1 — 1) (1 — e^^"^) sin (ttj/)] . We take the forcing term / corresponding to the exact solution. 



m 

in 



Figure 5.3(b) presents the streamlines of the approximate solution obtaine d by using the Argyris 
lent on a mesh with h — and 28550 DoFs. We note that Figure 5.3(b) resembles Figure 10 in 
Table 5.5 presents the errors cq, ei, and 62 for various meshsizes h. We note that the errors 



Table 5.5 follow the theoretical rates of convergence predicted by the estimates (4.21 )-(4.23 ) in 



Theorem 4.4 Again, we see that the orders of convergence in Table 5.5 are close to the theoretical 
ones for the fine meshes, but not as close for the coarse meshes. As stated previously, we attribute 
this to the inaccuracies at the thin (western) boundary layer at x = 0. 

5.2.3. Quasi-Geostrophic Equations. This section presents results for the FE discretiza- 



tion of the streamfunction formulation of the QGE (2.6) by using the Argyris element. To solve 



the resulting nonlinear system of equations, we use Newton's method with the following stopping 
criteria: the maximum residual norm is 10~*, the maximum streamfunction iteration increment is 
10"**, and the maximum number of iterations is 10. Our computational domain is = [0, 3] x [0, 1], 
the Reynolds number is Re — 1.667, and the Rossby number is Ro — 10~^. For completeness, we 
present results for two numerical tests, denoted by Test 5 and Test 6, corresponding to the exact 
solutions given in Test 1 and Test 2 of |3], respectively. 

Test 5: In this test, we take the same exact solution as that in Test 1 of [4], i.e., tjj{x,y) = 
sin^ (Try). The forcing term and homogeneous boundary conditions correspond to the 
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sm 



V 3 / 



Fig. 5.2. Linear Stommel Model \5.2\ , Test 2 f30f: Streamlines of the approximation, ip'^ , on a mesh with 
h = J- and 9670 DoFs. 



h 


DoFs 


eo 


L2 order 


ei 


ffi order 


62 


order 


1/2 


170 


0.00299 




0.04084 




0.7624 




1/4 


550 


3.217 X 10"^ 


6.539 


0.001031 


5.308 


0.04078 


4.225 


1/8 


1958 


3.437 X 10"^ 


6.548 


2.491 X lO"'' 


5.371 


0.002253 


4.178 


l/l6 


7366 


4.571 X 10"'' 


6.232 


7.026 X 10"'' 


5.148 


0.0001344 


4.067 


1/32 


28550 


6.704 X 10"" 


6.091 


2.113 X 10"** 


5.056 


8.26 X 10"" 


4.024 



Table 5.4 



Linear Stommel-Munk Model l |5.1[ l, Test 3 J^; The errors eo, e\, 62 for various meshsizes h. 




0.5 1 1.5 2 2.5 30 0.5 1 1.5 2 2.5 3 



(a) Test 3 [1. (b) Test 4 g|. 

Fig. 5.3. Linear Stommel-Munk Model ( |5.1| .• Streamlines of the approximation, i/i^ , on a mesh with h = ^ 
and 28550 DoFs. 



h 


DoFs 


eo 


L2 order 


ei 


order 


62 


H'' order 


1/2 


170 


0.06036 




1.162 




38.99 




1/4 


550 


0.01132 


2.414 


0.3995 


1.541 


21.4 


0.8656 


1/8 


1958 


0.0008399 


3.753 


0.05914 


2.756 


5.656 


1.92 


V16 


7366 


2.817 X 10"^ 


4.898 


0.004008 


3.883 


0.7378 


2.939 


1/32 


28550 


5.587 X 10"^ 


5.656 


0.0001607 


4.641 


0.0597 


3.627 



Table 5.5 

Linear Stommel-Munk Model ( |5.1[ l, Test 4 \4^: The errors eg, ei, 62 for various meshsizes h. 
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h 


DoFs 


eo 


L2 order 


ei 


ffi order 


62 


H'' order 


1/2 


170 


0.005709 




0.06033 




1.087 




1/4 


550 


3.726 X 10"^ 


7.259 


0.001086 


5.796 


0.04113 


4.724 


1/8 


1958 


3.597 X 10"^ 


6.695 


2.534 X lO"'^ 


5.421 


0.002252 


4.191 


1/16 


7366 


4.648 X 10"'' 


6.274 


7.065 X 10"'' 


5.165 


0.0001344 


4.067 


1/32 


28550 


6.737 X 10"" 


6.108 


2.116 X 10"** 


5.061 


8.26 X 10"" 


4.024 



Table 5.6 

QGE \2.G\, Test 5: The errors cq, ei, 62 for various meshsizes h. 




0.5 1 1.5 2 2.5 3 0.5 1 15 2 2.5 3 



(a) Test 5. (b) Test 6. 

Fig. 5.4. QGE \2.&\ : Streamlines of the approximation, , on a mesh with h = 7^ and 28550 DoFs. 



exact solution. 

Figure 5.4(a)| presents the streamlines of the approximate solution obtained by using the Argyris 



element on a mesh with h—^ and 28550 DoFs. We note that Figure 5.4(a) resembles Figure 7 in 



[30j . Table 5.6] presents the errors eo, ei, and 62 for various meshsizes h. The errors in Table 5.6 
follow the theoretical rates of convergence predicted by the estimates ( 4.21[ )-(4.23) in Theorem 4.4 



Again, since the exact solution does not display any boundary layers, we see that the orders of 
convergence in [Table 5.6| are close to the theoretical ones for the fine meshes, but are higher than 
expected for the coarse meshes. 

Test 6: In this test, we take the same exact solution as that in Test 2 of 2j, i.e., ipi^iU) = 
[(1 — 1) (1 — e~^°^) sin(7ry)]^. The forcing term and the homogeneous boundary conditions cor- 
respond to the exact solution. 



Figure 5.4(b) presents the streamlines of the approximate solution o btained by using the Argyris 
element on a mesh with h = ^ and 28550 DoFs. We note that Figure 5.4(b) resembles Figure 10 
in [30] . Table 5.7 presents the errors cq, ei 



and 62 for various meshsizes h. The errors in Table 5.7 
follow the theoretical rates of convergence predicted by the estimates ( 4.21[ )-(4.23) in Theorem 4.4 



We see that the orders of convergence in Table 5.7| are close to the theoretical ones for the fine 
meshes, but not as close for the coarse meshes. We attribute this to the inaccuracies at the thin 
boundary layer at x = 0. 

6. Conclusions. This paper introduced a conforming FE discretization of the streamfunc- 
tion formulation of the stationary one-layer QGE based on the Argyris element. For this FE 
discretization, we proved optimal error estimates in the i?^, and norms. A careful numerical 
investigation of the FE discretization was also performed. To this end, the QGE as well as the linear 
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h 


DoFs 


eo 


1/2 order 


ei 


order 


62 


H'' order 


1/2 


170 


0.3497 




1.9 




44.05 




1/4 


550 


0.0302 


3.533 


0.4279 


2.15 


21.74 


1.019 


1/8 


1958 


0.001507 


4.324 


0.06085 


2.814 


5.661 


1.941 


l/l6 


7366 


3.225 X lO"'^ 


5.547 


0.004042 


3.912 


0.7379 


2.94 


1/32 


28550 


5.672 X lO"'^ 


5.829 


0.000161 


4.65 


0.0597 


3.628 



Table 5.7 



QGE \2.6\, Test 6: The errors eo, ei, 62 for various meshsizes h. 



Stommel and Stommel-Munk models (two standard simplified settings used in the geophysical fluid 
dynamics literature 1301 [55] ) were used in the numerical tests. Based on the numerical results 
from the six tests considered, we drew the following two conclusions: (i) our numerical results are 
close to those used in the published literature [U [33 [33; and (ii) the convergence rates of the 



numerical approximations do indeed follow the theoretical error estimates in Theorems 4.1 and 4.4 



The convergence rates followed exactly the theoretical ones in the test problems where the exact 
solution did not display a thin boundary layer, but where somewhat lower than expected in those 
tests that displayed a thin western boundary layer, as expected. 

We plan to extend this study in several directions, including the time-dependent QGE and the 
two-layer QGE. 
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